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■ ABSTRACT 
C^-l \ The evolution of strange mode instabilities into the non linear regime has been followed 

by numerical simulation for an envelope model of a massive star having solar chemical 
composition, M = 50Af©, T c ff = 10 4 K and L = 1.17 • 10 6 L Q . Contrary to previously 
iyQ • studied models, for these parameters shocks are captured in the H-ionisation zone and 

, perform rapid oscillations within the latter. A linear stability analysis is performed to 

C ' verify that this behaviour is physical. The origin of an instability discovered in this way 

is identified by construction of an analytical model. As a result, the stratification turns 
out to be essential for instability. The difference to common stratification instabilities, 



, e.g., convective instabilities, is discussed. 

■ Key words: hydrodynamics - instabilities - shock waves - stars: mass-loss - stars: 

oscillations - stars: variables: other. 
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i 1 INTRODUCTION 

• »"H , Mass i ve stars are known to suffer from strange mode instabilities with growth rates in the dynamical range jKiriakidis. Fricke fc Glatzell 
/\ ' 1199.4 iQlatzel fc Kiriakidi"slll993h . The boundary of the domain in the Hertzsprung-Russel diagram (HRD) where all stel- 
^ ■ lar models are unstable - irresp ective of their metallicity -, coincides with the observed Humphreys-Davidson (HD) limit 
ijHumphrevs fc Davidson"1ll979h . Moreover, the range of unstable models covers the stellar parameters for which the LBV 
(luminous blue variable) phenomenon is observed. 

The high growth rates of the instabilities indicate a connection to the observed mass loss of the corresponding objects. 
To veri fy this suggestion, simulations of their evolution in to the non linear regime have been performed. In fact, for selected 
models Iciatzel. Kiriakidis. Chernigovskii fc Frickel (^999) found the velocity amplitude to exceed the escape velocity (see, 
however. iDorfi fc Gautschv I (|2000)). 

In this paper we report on a stellar model, which in the HRD is located well above the HD -limit, however, at lower 
effective temperature than the model studied bv loiatzel. Kiriakidis. Chernigovskii fc Frickel l)l999h . As expected, this model 
turns out to be linearly unstable with dynamical growth rates. When following the non linear ev olution of the i nstabilities, 
shocks form in the non linear regi me. The latter is customary in pulsating stellar enve lopes (see, e.g.. lChristvl lll966i) '). Contrary 
to the "hotter" model studied bv lGlatzel. Kiriakidis. Chernigovskii fc Frickel (Il999lh however, these shocks are captured by 
the H-ionisation zone after a few pulsation periods. The captured shock starts to oscillate rapidly with periods of the order 
of the sound travel time across the H-ionisation zone, while its mean position changes on the dynamical timescale of the 
primary, strange mode instability. This phenomenon is described in detail in section 3.1. Assumptions and methods on which 
the calculations are based are given in section 2. We emphasise, that in this publication, we concentrate on the oscillations of 
the captured shock. The phenomenon of shock capture by H-ionisation itself is not investigated here and will be studied in a 
separate paper. 

Apart from a detailed description of the shock oscillations found by numerical simulation the aim of the present paper 
consists of identifying their origin. This will be achieved by a linear stability analysis in section 3.2. It excludes a numerical 
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origin and attributes the oscillations to a secondary high frequencey instability in the shock zone. To identify the physical 
origin of the instability an analytical model is constructed in section 4. Our conclusions follow. 



2 BASIC ASSUMPTIONS AND METHODS 

2.1 Construction of initial model 

We investigate a stellar model having the mass M = 50Mq, chemical composition X = 0.7, Y — 0.28, Z — 0.02, effective tem- 
perature T fr = 10 4 A' and luminosity L — 1.17 ■ 10 6 Lq. These parameters have been chose n to ensure instability of the model . 
In the Hertzsprung- Russell diagram (HRD) it lies within the instability region identified bv lKiriakidis. Fricke fc Glatzell l|l993h 
(c.f. their figure 2). As only the envelope is affected by the instability, the model was constructed by standard envelope integra- 
tion using the parameters given above. The stellar core and nuclear energy generation are disregarded. Convection is treated in 
the standard mixing-length theory approach with 1.5 pressure scaleheights for the mixing length. T he onset of convection was 
determined by the Schwarzsch ild criterion. For the opacities, the latest versions of the OPAL tables ijlglesias. Rogers fc Wilsonl 
ll992llRogers fc Iglesiaslll992h have been used. 

2.2 Linear stability analysis 

Having constructed a hydrostatic envelope model its stability with respect to infinitesimal, spherical perturbations is tested. 
The relevant equatio ns corresponding to mass, en ergy and momentum conservation and the diffusion equation for energy 
transport are given in lBaker fc Kippenhahnl l)l962l) (hereafter BKA): 

(' = C* 4 (3C + C 5 p - C 6 t) (1) 

I' = (ia)C 1 {-p + C 2 t) (2) 

p' = _(4 + C 3 a 2 )C-p (3) 

t' = C 7 (^4C + C 13 l + C s p~C 9 t) (4) 

f, I, p and t are the relative perturbations of radius, luminosity, pressure and temperature, respectively, and dashes denote 
derivatives with respect to lnpo- cr is the eigenfrequency normalized to the inverse of the global free fall time rg = v/ R 3 /3GM. 
The coefficients d are determined by the background model where C13 denotes the ratio of total and radiative lumin osity. 
The o ther cofficients are defin ed in BKA. For the general theory of linear non-adiabatic stability, we refer the reader to ICoxl 
l|l980h and lUnno etafl Jl989h . 

The coupling between pulsation and convection is treated in the standard frozen in approximation, i.e., the Lagrangian 
perturbation of the convective flux is assumed to vanish. This is justified since the convective flux never exceeds 10% of the 
total energy flux. Moreover, the convective timescale is much longer than the dynamical time scale of the pulsations cons idered. 
The solution of the perturbation problem has been determined using the Riccati method dCautschv fc Glatzel lll99ol) . As a 
result of the linear non adiabatic (LNA) stability analysis we obtain periods and growth or damping rates of various modes 
together with the associated eigenfunctions. 

2.3 Non-linear evolution 

Having identified an instability by the LNA analysis its growth is followed into the non-linear regime. Assuming spherical 
symmetry, we adopt a Lagrangian description and choose as independent variables the time t and the mass m inside a sphere 
of radius r. The evolution of an instability is then governed by mass conservation, 

— -— = (5) 

dm 4-rrp 



momentum conservation, 
A 2 dP Gm 

+ ^ r \z: + ^r = (6) 



energy conservation, 
dL_Pdp 9_E 

dm p 2 dt dt { ' 

and the diffusion equation for energy transport, 

dT 3k(L — Lkonv) _ ~ /o\ 
~dm~ ~ 6ATY 2 acr 4 T 3 ~ [ ' 
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Table 1. Unstable modes of the initial model. 

o> 0.53 1.22 1.66 2.12 2.26 3.34 3.86 
<Ti -0.06 -0.18 -0.13 -0.20 -0.12 -0.04 -0.04 

a r denotes the real part and ai the imaginary part of the eigenfrequency <r normalized by the global free fall time. 

where p, p, T, L, and E denote density, pressure, temperature, luminosity and specific internal energy, respectively, a is the 
radiation constant, c the speed of light and G the gravitational constant. For consistency, the equation of state p(p, T) and the 
opacity k are identical with those used for the construction of the initial model. In accordance with the LNA stability analysis, 
convection is treated in the frozen in approximation, i.e., Lk onv is taken to be constant during the non-linear evolution and 
equal to the initial value. For the treatment of shocks artificial viscosity is introduced by substituting P = P + Q with (v is 
the velocity) 

(C p(divv) 2 divu<0 
\o div v > 

and Co > (von Neumann - Richtmyer form of artificial viscosity). 

For some difference schemes including the Fraley scheme, which th e present method is based on, this form of the artificial 
viscosity can give rise to undesi red, unphysical oscillations (se e, e.g., iBuchler fc Whaler] lll990t0 . To avoid these, artificial 
tensor viscosity is usually used l|Tscharnuter fc Winklerl jl97flh L In order to be sure that the oscillations observed are not 
caused by the form of the artificial viscosity, we have run tests both with volume and tensor viscosity. As a result, shock 
oscialltions are found independently for any form of the artificial viscosity. As the von Neumann - Richtmeyr viscosity allows 
for a straightforward formulation of the boundary conditions discussed below, we have for convenience chosen to work with 
it. 

The inert hydrostatic core provides boundary conditions at the bottom of the envelope by prescribing its time independent 
radius and luminosity there. As the outer boundary of the model does not correspond to the physical boundary of the star, 
boundary conditions are ambigous there. We require the gradient of heat sources to vanish there: 

grad(divf) = (10) 

This boundary condition is chosen to ensure that outgoing shocks pass through the boundary without reflect ion. The nume rical 
code relies on a Lag r angia n, with respect to time implicit, fully conservative difference scheme proposed bvlFral ev I Jl96 Sl) and 
ISamarskii fc Popov! dl969l) Concerning tests of the code, we adopted the same criteria as lClatzel. Kiriakidis. Chernigovskii fc Frickel 
Il999h . 



3 NUMERICAL RESULTS 

3.1 The evolution of the stellar model 

Density and temperature of the initial model as a function of relative radius, are shown in figures 0al-0a2. The stratification 
exhibits a pronounced core-envelope structure, which is typical for stellar models in this domain of the HRD. More than 96 
per cent of the mass is concentrated in the core, which extends to less than 5 per cent of the total radius. It remains in 
hydrostatic equilibrium and is not affected by the instability. 

The initial model has been tested for stability and been found unstable on dynamical timescales. This instability will 
be referred to as primary instability hereafter. With respect to its physical origin, it is a strange mode instability, which has 
been identified in a variety of stars including Wolf-Rayet-stars, HdC-stars and massive stars (like the present model). Strange 
modes appear as mode coupling phenomena with associated instabilities whenever radiation pressure is dominant. The latter 
is true for a large fraction of the radius in the present model. The linear stability analysis of the initial model reveals several 
unstable modes. Eigenfrequencies of the most unstable ones, i.e., their real (oy) and imaginary parts (cr;), are presented in 
table □ 

The evolution of the linear instabilities was followed into the non-linear regime by numerical simulation using the hydro- 
static model as initial condition. No additional initial perturbation of the hydrostatic model was added. Rather the code was 
required to pick the correct unstable modes from numerical noise. By comparing growth rates and periods obtained in the 
simulation with the results of the LNA analysis, the linear regime of the evolution was used as a test for the quality of the 
simulation. 

In the non-linear regime sound waves travelling outwards form shocks and initially inflate the envelope to 2.5 initial radii. 
Thus velocity amplitudes of 10 7 [cm/sec] are reached. One of the subsequent shocks is captured around the H- ionization zone 
at relative radius r/R = 0.58 and 3.6 < logT < 4.7. The mechanism responsible for the shock capturing will not be studied 
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Figure 1. Temperature T and density p as a function of relative radius of the initial model (al-a2), and the model at 5 ■ 10 7 sec (bl-b2). 
Velocity v and Mach number M are shown as a function of relative radius for the model at 5 ■ 10 7 sec in (c) and (d), respectively 



in this publication. Rather, we will investigate the oscillations on the shock front and show that they are of physical origin. 
A snapshot at time t = 5 ■ 10 7 sec of the situation containing the captured shock is shown in Figures^bl and0b2. Figure^c 
shows the velocity as a function of relative radius at this instant. Sound waves are generated in the region around r/R m 0.1 
and travel outwards, growing in amplitude and steepening. In the snapshot one such wave is located at r/R ss 0.25. The 
captured shock front is located at r/R ~ 0.58 and the outer envelope is collapsing onto it. The small panel in Figure Qc 
shows the details of the region containing the captured shock, indicating the grid resolution by (x). Within the Lagrangian 
description, ~ 150 of the 512 gridpoints used are concentrated in the shock zone. Figure^d shows the Mach number M — v/v 3 
as a function of relative radius for the snapshot (v s is the local sound speed). The Mach number changes by 3.5 across the 
shock front around r/R~ 0.58. 

After the formation of the captured shock its position varies only weakly by w 0.2 relative radii on the timescale of the 
primary instability (see figure 0b). Superimposed on this variation is a much faster oscillation, whose timescale is related to 
the sound travel time across the shock (w 10 5 sec). It is even more pronounced in the run of the luminosity (figure |3Jc2). 
The onset of the fast oscillation with the capturing of the shock by the H-ionization-zone is illustrated in figures 0a andEJcl, 
where the velocity v and relative luminosity L/Lo at the outer boundary are shown as a function of time. Up to « 4 • 10 7 sec 
the velocity varies on the timescale of the primary instability and the luminosity remains approximately constant due to the 
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Figure 2. The velocity v at the boundary (a), the relative position of the shock front (b) and the relative luminosity L/Lq at the 
boundary on two different scales (cl-c2) as a function of time t. 



low heat capacity of the envelope of the star. After « 4 • 10 7 sec, when the shock has been captured by the H-ionization-zone, 
luminosity and velocity vary on the shorter timescale of the secondary shock oscillation. The luminosity perturbation has its 
origin in the shock. Due to the low heat capacity the luminosity perturbation remains spatially constant above the shock. 

In principle, the high-frequency secondary oscillations of the shock could be caused numerically. However, the results are 
largely independent of the numerical treatment and parameters, which has been veryfied by extensive numerical experiments 
suggesting a physical origin of the phenomenon. In section 13.21 we shall argue in favour of the latter by presenting a linear 
stability analysis providing an instability with appropriate frequencies and growth rates. 



3.2 Stability analysis of a model containing a captured shock 

In this Section we shall initially assume and then prove a posteriori, that the secondary oscillations of the captured shock 
described in section 13.11 are caused by physical processes. We perform a linear stability analysis of a background model 
by assuming that the dependent variables radius, pressure, temperature and luminosity may be expressed as the sum of a 
background contribution and a small perturbation: 

x(m, t) = xo(m, t) + Xi(m, t) for x £ {r, p, T,L} (11) 

The background coefficients xo(m, t) may be regarded as time independent, i.e., xo(m, t) = xo(m), as long as the perturbations 
vary on much shorter timescales as the background, i.e., as long as the condition 

rflogso(m,t) dlogxi(m,t) , . 

dt * dt ( ' 

holds. 4? denotes the Lagrangian time derivative. Thus, the variations on dynamical timescales of the model containing the 
captured shock are regarded as stationary with respect to the anticipated much faster instability. Eigenmodes with periods 
of the order of the dynamical timescale suffer from the competition with the variation of the "background" model, whereas 
the approximation holds for those with much shorter periods. For the model considered, the approximation is correct for 
\a\ > 100. Should unstable eigenmodes of this kind exist, this would prove the instability and the high frequency oscillations 
of the captured shocks to be of physical origin. Therefore, the results of a linear stability analysis of such a model are 
meaningful, as long as the obtained frequencies are interpreted properly. 
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A problem with this strategy is that the numerical simulations provide only the superposition of the slow dynamical 
and the secondary, fast oscillations. The linear stability analysis, however, requires a - on the fast oscillations - stationary 
background model. It is obtained by an appropriate time average over a numerically determined sequence of models. "Ap- 
propriate" means, that the average has to be taken over times longer than the short period oscillations and shorter than the 
dynamical timescale. Thus all physical quantities Q(m,t) are averaged according to 

<Q(m)>= _ / Q(m,t)dt (13) 
t e t s Jt 

where t s and t e are the beginning and the end of the averaging intervall and satisfy the requirements discussed above. t s has 
been varied between 4 • 10 7 sec and 5 • 10 7 sec (after the formation of the shock front) and the averaging interval between 5 • 10 5 
sec and 1 • 10 6 sec. All averages exhibit qualitatively the same behaviour and the LNA stability analysis is largely independent 
of the averaging parameters. The results presented in the following were obtained for t s = 5 • 10 7 sec and t e = 5.05 • 10 7 sec. 

With these assumptions, the linear perturbation equations 11141 which have been derived for a strictly static background 
model, remain valid even for the situation studied here, except for the momentum equation [3] which has to be modified 
according to: 

p> = -(4 + C 3 a 2 )(~Q lP with Q 1 = _|gigi (14) 

Qi 7^ 1 accounts for the deviations from hydrostatic equilibrium. 

As a result of the linear stability analysis (with Q\ 7^ 1) , the expected unstable modes having high frequencies have been 
identified. E.g., a typical mode of this kind satisfying the assumptions discussed has the frequency oy = 4837.2 and the growth 
rate <Ji = —83.9. 

In a second step, we investigate the influence of deviations from hydrostatic equilibrium, i.e., the deviations from Qi = 1. 
For this purpose we rewrite equation 1141 as 

p' = -(4 + C 3 a 2 )C-p + ^(l-Q 1 )p (15) 

with ^ <f> ^ 1. The limits $ = and $ = 1 correspond to hydrostatic equilibrium and the averaged model containing 
the shock, respectively. The influence of deviations from hydrostatic equilibrium is then studied by varying $ between and 
1. Following the mode having ay = 4837.2 and a\ = —83.9 at $ = 1 to $ = its frequency and growth rate changes to 
o T = 4651.7 and a\ = —119.6. The moduli of the corresponding Lagrangian displacements, which indicate the kinetic energy 
of the pulsations, are shown in Figure as a function of relative radius. The energy of the pulsation is concentrated around 
r/Rm 0.58 and drops off exponentially above and below. As a result, neither eigenvalues nor eigenfunctions differ significantly 
for $ = and $ = 1, i.e., the assumption of hydrostatic equilibrium is justified for the unstable modes considered. Therefore, 
we will assume hydrostatic equilibrium in a further discussion and investigation of the secondary instability, i.e., all subsequent 
results were obtained assuming $ = 0. 

The results of a LNA stability analysis according to eouations lll4l and Section ^. 21 for the averaged model are summarized 
in Table|5] where representative values for the eigenfrequencies of unstable modes are given. Three sets of unstable modes may 
be distinguished. Low order modes with ay between 0.9 and 9 have growth rates of the order of 0.2, i.e a ratio of ^- w 0.1. 
They can be identified with the primary instability. However, their periods compete with the variation of the background 
model and therefore these modes have to be interpreted with caution. The properties of two classes of high order unstable 
modes with frequencies ay between 140 and 4650 are in accordance with our approximation. One of them has high growth 
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Table 2. Unstable modes of the averaged model 
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-26.3 
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-0.35 
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ay denotes the real part and ai the imaginary part of the eigenfrequency <r normalized by the global free fall time. 

rates with a ratio of ^ ~ 0.1, the second low growth rates with a ratio of ^p- ~ 5 • 10~ 4 . The latter may be identified with high 
order primary instabilities, whereas the former are attractive candidates for the secondary, shock front instabilities sought. 

For further discussion, we consider eigenf unctions and the corresponding work integrals of representative members of the 
different s ets of modes. The work integral is a widely used tool to identify the regions in a star, which drive or damp the 
pulsation, loiatzell ljl994h has shown, that the concept of the work integral is not necessarily restricted to small values of the 
damping or growth rate. By replacing the conventional time average by an ensemble average it can be extended to arbitrary 
values of In any case, one arrives at the expression 



for the work integral, where Im(z) denotes the imaginary part of z, ()* denotes complex conjugation and p\ denote the 
spatial parts of the eigenfunctions of the relative pressure and density perturbations, respectively, p is the pressure of the 
background model. The sign of the integrand in eg uat ion 1 1 6 1 det ermines . if a region of the star damps or drives the pulsation, 
where lm(p^ p^*) < corresponds to driving and Im(p^ p^*) > to damping. Some authors (e.g. BKA) use logp instead of r as 
independent variable and therefore obtain an opposite sign of the differential work integral for driving and damping influence. 
To match this convention, —W(r) is shown in figures^Jb, i.e., driving regions correspond to positive —W(r), damping regions 
to negative — W(r). 

According to its eigenvalues the mode corresponding to ay = 2.22, ay = —0.33 was identified as a primary instability. 
This is supported by the Lagrangian displacement component £ of the eigenfunction and the work integral shown in figure 
2]al and^Jbl. The shock front acts as an acoustic barrier causing the eigenfunction to vanish above it (figure 2]al). The 
work integral (figure0]bl) exhibits two driving regions which coincide with the opacity peaks at logT = 5.3 (caused by the 
contributions of heavy elements) and logT = 4.7 (He-ionization). The stability properties are not affected significantly by the 
shock region. 

The counterparts of figures 0]al and0]bl for a weakly unstable high frequency mode having a r = 157.2, Oi = —0.08 are 
shown in figures Ha2 and|lb2. A gain, the isolating effect of the shock front causes the dramatic variation of the amplitude 
around r/Ro = 0.6. However, in contrast to the eigenfunction presented in figure^Jal the amplitude is now significant above 
and negligible below the shock. High order modes of this kind in general exhibit strong damping. For the mode considered 
the shock efficiently screens the inner damping part of the stellar envelope. Thus the region below the shock contributes only 
weak damping which is overcome by the driving influence of the shock as shown by the work integral in figure |l]b2. 

Apart from splitting the acoustic spectrum by an acoustic barrier into two sets of modes associated with the acoustic 
cavities below and above the shock, respectively, the shock itself gives rise to a third set. Lagrangian displacement and work 
integral for a typical member of this set having a r = 162.9, at = —35.4 are shown in figures 0]a3 and0|b3. The amplitude 
of this unstable mode reaches its maximum on the shock and drops off exponentially above and below. Note its oscillatory 
behaviour on and close confinement to the shock. The real parts of this set of eigenfrequencies of ~ 200 — 500 correspond to 
periods of n w 8 • 10 4 - 2 ■ 10 5 sec, which are observed in the luminosity perturbations (cf. figure |5]c2) induced by the shock 
oscillations. The work integral (figure 0]b3) shows, that the shock is driving this instability, and that the regions above and 
below do not contribute. Moreover, the basic assumption of stationarity of the averaged model holds for the frequencies and 
growth rates obtained. 

Thus we have identified an instability by linear analysis of an averaged model, which resembles the shock oscillations 
observed in the numerical simulations, both with respect to timescales and spatial structure. We therefore conclude, that 
the shock oscillations are not numerical artifacts. Rather they have a physical origin and are caused by an instability whose 
mechanism will be investigated in detail in the following sections. 

3.3 Approximations 




(16) 



In order to gain further insight into the physical processes responsible for the instability, different approximations in equations 
have been considered. To obtain a continuous transition from the exact treatment to the approximation, we introduce 
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Figure 4. Lagrangian displacements f (a) and integrated workintegrals (b) as a function of relative radius for the eigenfrequencies 
oy = 2.22, Ui = -0.33 (1), o> = 157.2, <jj = -0.08 (2) and oy = 162.9, Ui = -35.4 (3) of the averaged model. 



a parameter $ with $ = 1 corresponding to the exact problem and <E> — > 0, oo to the approximation. The numerical results, 
i.e. the eigenvalues of the shock front instabilities, are followed as a function of $. 
Introducing <E> into the Euler equation as 

p' = -(4 + $-C 3 a 2 )C-p (17) 

the limit <3> — ► corresponds to vanishing acceleration and implies the elimination of acoustic modes from the spectrum, 
which then only consists of secular modes. Application of this limit to the shock instabilities has not revealed any unstable 
modes. Rather the eigenvalues have diverged. This excludes a thermal origin of both the unstable modes and the instability 
mechanism. For a proper treatment of the instability, the mechanical acceleration has to be taken into account. 
Introducing <!> into the equation for energy conservation as 

I' = C 1 -Q-{ia)(-p + C a t) (18) 



the adiabatic limit is obtained by $ — > oo. The latter implies (— p + C-zt) = 0, i.e. the algebraic adiabatic relation between 
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Figure 5. The coefficients C5 = a (a) and C7 = V (b) of the averaged model as a function of relative radius. 



pressure and temperature perturbation. No unstable modes have been found following the shock instabilities into the adiabatic 
limit. 

Introducing <3> into the equation for energy conservation as 

I' = Ci ■ $• (ia)(-p + C 2 t) (19) 

$ —* corresponds to the so called NAR-limit (Non- Adiabatic- Reversible limit) (|Gautschv fc Glatzel lll99(l l. Although this 
approximation - like the adiabatic approximation - implies constant entropy, it does not represent the adiabatic limit (— p + 
GiC) = 0. Rather it is equivalent to Ci — > 0. Since Ci is related to the thermal and dynamical timescales r t h and Td yn by 



Tth 
7~dyn 



Vac 



C3C4 



(20) 



this approximation is also being refered to as the zero thermal timescale approximation (r and V a d are the adiabatic indices) . 
Physically, it means that the specific heat of the envelope is negligible and luminosity perturbations cannot be sustained. In 
particular, this approximation rules out the classical jc-mechanism as the source of an instability - should it exist in the NAR- 
limit - since this Carnot-type process relies on a finite heat capacity. When following the frequencies of the modes belonging 
to the shock front instabilities into the NAR-limit, periods and growth rates change only slightly (by at most 10 per cent). 
Thus the NAR-approximation may be regarded as a satisfactory approximation and will form the basis of our investigations 
in the following sections. 



4 AN ANALYTICAL MODEL 
4.1 Three-Zone-Model 

The modal structure identified in section with three sets of modes associated with three acoustic cavities (inner envelope, 
shock and outer envelope) suggests the construction of a three zone model. In order to enable an analytical solution, the 
coefficients of the differential equations are kept constant in each zone. 

According to the previous section the NAR-approximatin is sufficient to describe the shock front instabilities. The equation 
of energy conservation is then satisfied identically and luminosity perturbations vanish. Thus we are left with a system of 
third order comprising the mechanical equations and the diffusion equation with zero luminosity perturbation. 

Further reduction of the order of the differential system is achieved by considering its coefficients which depend on the 
properties of the averaged model. In figure the coefficients C5 = a = §j§|^ | ~ ^ and O7 = V = diogp are snown as a 

function of relative radius. f3 denotes the ratio of gas pressure to total pressure. The coefficients C4 = — j|° sr and C3 = 4 ^ ? 
may be regarded as constant all over the envelope. Approximate values are C4 ~ | and Cs ~ 3. The latter holds because 
almost the entire mass is concentrated in the stellar core. From figurc|^]we deduce that radiation pressure is dominant except 
for the shock zone. Therefore we replace the diffusion equation^Jby the algebraic equation of state for pure radiation (p = At) 
in the inner and outer envelope. On the other hand CV = V can - to first approximation - be regarded as singular in the shock 
zone. According to equation ^ this requires the expression (— 4£ + Csp — Cgt) to vanish there. Thus the differential diffusion 
equation is replaced by an algebraic re lation in all three zones, red ucing the system to second order. 

Adopting the alternative notation JlBaker fc Kippenhahnll962l) Ce = 5, Cs — n P and C9 = 4 — kt, where 5 is the negative 
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logarithmic derivative of density with respect to temperature at constant pressure, k p the logarithmic derivative of opacity 
with respect to pressure at constant temperature and kt the logarithmic derivative of opacity with respect to temperature 
at constant pressure, and choosing the relative radius x as the independent variable, we are left with the following set of 
equations: 

Yi = i(3C + ap-Jt) (21) 
\t = -(4 + 3* 2 K-, (22) 




t = <v ,-, Ty («PP-4C) (23) 

x G [0, a) or x G (b, 1] 

a and b denote the lower and upper boundary of the shock zone. The transformation of the independent variables lnpo - *• x 
introduces the factor tp, which is constant within the framework of the three-zone-model, and given by an appropriate mean 
of the quantity . In general ip is negative and of order unity. 

We are thus left with a system of second order consisting of the mechanical equations (continuity and Euler equations) 
which is closed by the algebraic relations = — 4£ + Cgp — C%t and p — 4 • t for the shock region and the inner/outer regions, 
respectively. We rewrite it as: 

dc 

dx 
dp 
dx 

where 



Ai, 2 • C + B 1>2 ■ p (24) 
C ■ C + D ■ p (25) 




x G [a, 61 

A h2 = { r V 4 -"W 1 ' J (26) 

x G [0, a) or x G (6, 1] 

(27) 

Lh , - { T V" 4 ~ K7 7 J (28) 

^(f-n) ie[o,a)one(M] 

(29) 

C =-V(4 + 3a 2 ) (30) 
D =-tp (31) 

and the subscript 1 denotes the values of the coefficients in the shock region, the subscript 2 values in the inner and outer 
regions. We introduce new variables by 

C = e J M - 2dx ■ C (32) 
p =J Ddx -p (33) 

The system I24l25l then reads 

^ 7-» ~ [ Ddx — f Ai ?dx /o *\ 

-r 1 = Sio-p-e J -e 1 1,2 (34) 

# = C ■ C ■ e~ 1 Ddx ■ e 1 Ai - 2dx (35) 
dx 

These equations are equivalent to the following single second order equation: 

d_ I J_ / Ddx — J A 12 dx dp 
dx \C dx 



d ( 1 ! Ddx - f Ai odx dp\ f odx - f Ai ->dx D n 

— 1 — • e J • e J 1,2 ■ — — e J e J 1,2 • Bi^ • p = (36) 



4-1.1 Mathematical Structure of the Problem 
Equation 1361 may be written as 

d ( f Dda: — f Ai adx \ , ^ / f Ddx - f Ai odx 7-, ~ . 2 / f Ddx - f Ai odx t-> * rt /ot\ 

— le J -e 1 1,2 ■ — I + 4 • i/j • e J e J 1,2 -iJi^-p + Scr • ip ■ e' e ' 1,2 • -Bi,2 ■ p = (37) 



and has the form 
d 



dx 



w(x)p + Xu{x)p = (38) 
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with q{x) > in the integration intervall. However, u(x) is positive in the inner and outer regions and negative in the shock 
region, i.e., u(x) changes sign in the integration intervall. (This holds also for w(x).) Therefore, this problem is not of Sturm- 
Liouville type. On the other hand, if we consider each zone separately with boundary conditions p — 0, equation I38l describes 
a Sturm-Liouville problem. In the shock zone we define eigenvalues A = —a 2 and thus have u(x) > 0, w(x) > 0, for the inner 
and outer zones we get u(x) > 0, w(x) < by defining A = a 2 . 

For a Sturm-Liouville problem, the eigenvalues are real and form a sequence 

Ai < A 2 < A 3 < A 4 < . . . (39) 

Furthermore, Ai may be estimated on the basis of the variational principle 

, . Iol{x)\p'\ 2 dx + J^w(x)\pfdx ,, n , 

Xl = m - m r i : ,,.,2 , ( 4 °) 

V j u{x) \p\ dx 

Therefore Ai = —o\ is positive in the shock zone, since w(x) is positive there. This means that we have purely imaginary 
eigenfrequencies Oj = ii-v/A^ with positive \j and 

\cri\ < \a 2 \ < \a 3 \ < | ff4 1 < • • • (41) 

Thus the shock zone provides unstable eigenfrequencies. 

Since w(x) is negative in the inner and outer zones, we cannot guarantee Ai to be positive there. For sufficiently large j, 
however, Xj will always become positive. As a consequence, all eigenfrequencies ttj = iA^Aj will become real for sufficiently 
high order j ^ n and satisfy: 

\(T n \ < |c«+l| < cr, x +2 1 < jcr„+3| < . . . (42) 

In principle, the mathematical structure of the problem allows for imaginary pairs of eigenfrequencies at low orders in the 
inner and outer zones. For the particular parameters studied in the following sections, however, Ai turned out to be positive, 
i.e., n = 1 and all eigenfrequencies are real. 

Even if equation 1371 together with the boundary conditions p = at x — and x = 1 (three-zone-model) is not of 
Sturm-Liouville type, the differential operator 

v =4:(<i4:) +*u-w (43) 



dx \ dx J 

in equation 1381 can be shown to be self adjoint with the boundary conditions p — at x = and x = 1. Therefore the 
eigenvalues A are real and we do expect only real or purely imaginary eigenfrequencies a, i.e., we will not be able to reproduce 
the complex eigenfrequencies of the exact problem in this approximation. 



4.1.2 Results 

Assuming the coefficients C and -E?i,2 to be constant, equations 1341 and 1351 are solved by the Ansatz p, £ oc e . For the 
wavenumbers k we get 



k = ±y/B ll2 -C (44) 
Thus the general solutions reads 

(^■e^Vx + az-e-^ 7 * x G [0, a) 
p = < h ■ e^ 1 + b 2 ■ e~' /BTUx x G [a, b] (45) 
[ Cl ■ e^ Ux + c 2 ■ e~ vl ^ Ux xe(b,l] 

a.1,2, 61,2 and ci,2 are integration constants and have to be determined by the requirements of continuity and differentiability 
of p at x = a and x = b and the boundary conditions at x = and x = 1. For the latter we choose p = 0, which implies 

<22 = —ai and (46) 

c 2 =- Cl -e 2 ^ (47) 

Together with the requirements of continuity and differentiability this yields the dispersion relation 

{vwc - ^mc)e^ Uib - i) + (-^mc - ^mc)e-^ U{b - i) 



(48) 
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Table 3. Eigenfrequencies a [a r : real part, cr;: imaginary part) of the three-zone-model having the parameters Bi = — 4i/>, B2 = ^, 
a = 0.57, b = 0.59, ip = -1 



<7 r 


12.01 


16.78 


18.55 


24.86 


25.83 


31.29 


34.67 


Ti 























<7 r 























7.21 


52.43 


97.77 


143.11 


188.46 







where the eigenfrequencies a are contained in the coefficient C. In general, the roots of equation 1481 have to be calculated 
numerically, using, for example, a complex secant method. Separate spectra for the three isolated zones may be obtained by 
assuming the boundary conditions p — at x — a, b instead of continuity and differentiability requirements. We are then left 
with the dispersion relations 

D = (VBic - vmc)e^ ja + (-vmc - vmc) e -^ Ua 
D = (vmn + vmc)e^ u ^ + (- jmc + VBic) e -^ U( - b - 1 ^ 

° Uter (V5r? - VftC)^*'- 11 + - ( ' 

D 8hock = 1 = e 2 vW(-*) (51) 

for the inner, outer and shock zones, respectively. 

For the averaged model we have Bi ~ —4ip and dominant radiation pressure implies B2 w \ip. Inserting these values into 
equations !49l5Ul we are left with 



)/£(4 + 3o> = A^+i2ZE nGZ (52 ) 



^(4 + ^(6-1) =^±^L 



neZ (53) 



Thus we have real a, i.e., neutrally stable modes, if the inner and outer regions are considered separately, in accordance with 
the discussion in section Pi. 1.1 1 For the shock region equation 1511 yields 



ai 




7x +a 2 


■ 




x G [0, a) 




61 




7x + b 2 






x g [a, b] 


(55) 


Cl 


■ 


7x + c 2 






x G (6,1] 





2^4V> 2 (4 + 3a 2 )(6 - a) = 2nni neZ (54) 

These solutions correspond to purely imaginary a implying instability. The solutions of equations 1521541 can be used as 
initial guesses for the numerical iteration of equation 1481 the disperion relation of the three-zone-model. Some representative 
eigenvalues of the three-zone-model are given in Table |31 

Once the eigenfrequencies are determined, the corresponding eigenfunctions are given by 

p = < 

The factor e^ x is due to the transformation from p to p. 

Typical eigenfunctions are presented in figures HJal-E]a3. Three types of modes may be distinguished belonging to the 
three zones of the model. Real eigenfrequencies are associated with the inner and outer region. Except for the shock region 
they are oscillatory and reach their maximum in the respective region. "Shock modes" correspond to unstable and damped 
modes (purely imaginary pairs of eigenvalues). They oscillate in the shock region and are evanescent elsewhere. We note the 
correspondence of figuresHJal andPj]al,|S|a2 andPj|a2 and|S|a3 andPj]a3, i.e., the results of the analytical model resemble those 
of the exact analysis. The influence of the shock position on the modal structure may also be studied within the framework 
of the three-zone-model. As long as the width of the shock zone and the coefficient B\ are not varied, the "shock modes" are 
not affected. The dependence on the shock position of the neutrally stable "inner" and "outer" modes is shown in figure|S|b. 
Moving the shock position outwards, the frequencies of the inner modes decrease, whereas those of the outer modes increase, 
according to the variation of the length of the corresponding acoustic cavities. This leads inevitably to multiple crossings 
betwe en the frequencies of the inner and outer modes, which unfold into avoided crossings (see, e.g., lOautschv fc Glatzel I 
(1990)). Mode interaction by instability bands is excluded here according to the general discussion in section l4.1.1l 



4-1.3 Interpretation 

The three-zone-model reproduces the effects of the shock front regarding important aspects: The front acts as an acoustically 
isolating layer which separates the inner and outer part of the envelope. As a result, these parts provide largely independent 
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Figure 6. Eigenfunctions for the three-zone-model with the parameters B\ = — 4i/>, B2 = -j, a = 0.57, b = 0.59, V = — 1) an d the 
frequencies a r = 44.04, Oi = (al), <r r = 43.53, at = (a2), <r r = 0, Oi = 188.46 (a3). (b): Eigenfrequencies a r of neutrally stable modes 
as a function of the position a of the lower boundary of the shock region for fixed b — a = 0.02 and B± = — 4i/>, B2 = i , if> = —1. 



spectra. This may be illustrated by the variation of the position of the shock front. Apart from the expected spectra associated 
with the inner and outer envelope, an additional spectrum of modes is generated by the shock region itself. 

Comparing eigenfunctions of the averaged and the analytical model (figures 2]al-a3 and figures |S|al-a3), we find a 
strikingly similar behaviour. In particular, the confinement of the unstable shock modes is present in both cases. Due to 
constant coefficients, however, the analytical model reproduces neither decreasing amplitudes nor increasing spatial frequencies 
towards the stellar center. 

We have identified unstable modes in the shock zone of the analytical model. They resemble those of the shock instabilities 
of the averaged model, and are related to the sound travel time across the shock zone. Its radial extent is primarily responsible 
for their high frequencies. 

The analysis in section l4.1.1l has shown, that the sign of u(x) in equation is responsible for the instability in the shock 
region. This sign is determined by the term —r~, which is given by 



Bi 



a 
3 



(56) 



4 — kt 

Estimating the various terms in equation 1561 we find that the sign of k v determines the sign of A dependence on the sign 
of k p of the instability, however, is not recovered in the exact problem, which can be tested by replacing n p with — k p there. 
The exact problem is not affected by this substitution. Thus we conclude, that the analytical model does not provide correct 
results in this respect and needs to be refined to describe the instability properly. In order to investigate the origin of the 
instability, some of the simplifying assumptions of the analytical model need to be dropped. In this direction, a more realistic 
model of the shock zone will be presented in the following section. 



4.2 Shock-Zone-Model 

Our study of the three-zone-model in section im has shown, that inner, outer and shock zones may to good approximation be 
treated separately by assuming suitable boundary conditions, e.g., vanishing pressure at boundaries and interfaces. Moreover, 
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the instabilities of interest are not provided by the inner and outer zones. Therefore we restrict the following study to the shock 
zone by applying the boundary conditions p(a) = p(b) — 0. Within the framework of the analytical model the coefficients of 
the perturbation equations are taken to be constant with the values given in section |4~T1 

Contrary to section l4~Tl we will not replace the diffusion equation by an algebraic relation here, as this turned out to lead 
to erroneous results. However, we still adopt the NAR-approximation. The set of equations considered then reads: 

ldC 1 , 
ip dx 
1 dp 
ip dx 
ldt_ 
ip dx 
1 dl 
ip dx 

Written in matrix form this yields 



= V 



-m + ap-St) 
(4 + 3a 2 )C-p 

4C + K P p- (4- K T )t) 







(57) 
(58) 
(59) 
(60) 



LA. 

ip dx 



p 
t 

W 



/ i 

-(4 + 3a 2 ) 
-4V 







_ 6 
3 



-V(4-« r ) 




o 





p 
t 

W 



(61) 



The differential equation is solved by an exponential dependence oc e %kx of the dependent variables. Thus we arrive at the 
linear algebraic equation 



(62) 



This equation has a non trivial solution only if the determinant of the matrix vanishes, which provides a quartic equation for 
the wavenumber k. One of its roots is zero, the remaining three roots are determined by the following cubic equation: 



/ 


1 ik 

— 


a S 
3 3 


o \ 


(C) 




f°\ 




-(4 + 3a 2 ) 


-i - f o 





p 









-4V 


Vk p -V(4-kt)-| 





t 







V 








ik 

•0 / 


W 




w 



1 



V(4- k t ) 



+ 



+ 



(4- k t ) V(4 — k t ) 



(1_^(4 + 3a 2 )) + 



-l + f(4 + 3a 2 ) 



U-k,) yr^ + ^ + p 



= 



(63) 



do 



In the limit of large V they may be given in closed form: 



d 1 dl 



(64) 



= -V(4-kt) 



(65) 



The general solution to the perturbation problem consists of a superposition of four fundamental solutions associated with 
the four roots for the wavenumber, two of which are oscillatory (those associated with fci and kz). The dispersion relation is 
then derived by imposing four conditions. In addition to the boundary conditions p — at x = a, b, we require the two non 
oscillatory fundamental solutions not to contribute to the eigensolution. The latter is then only determined by fci and kz: 



7 Zfcl X I 7 ikn X 

p = hi ■ e + hi ■ e 



(66) 



where h\ and hi are integration constants. They are determined by the boundary conditions p — at x — a, b, which imply 
2irn 



ki — k 2 



(a-b) 

where n € Z denotes the order of the overtone. Using equation 1641 we get 

dl ~ 2 ~ 2 



(67) 



di 



ip 2 (a~b) 2 



(68) 



With the definitions of di and d? (equation I63H we arrive at a quadratic equation in a 2 . Expanding the coefficients of a 2 
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terms of -i- and assuming 



2 a 2 



V 2 (4-k t ) 2 
Defining 

- _ V 2 (4- KT ) 2 



WW 
5 k 



to be large, we obtain to lowest order in ^: 



+ 



(4-Kt) J tp 2 {a-b) 



2a 2 

6 



this equation has the solutions 



5k h 



(4- KT ; 



2 



(4-kt; 



V ^ 2 (a-b) 2 



(4-«r) 



(69) 



(70) 



(71) 



In the NAR-approximation, eigenfrequencies come in complex conjugate pairs, i.e., complex eigenfrequencies imply instability. 
According to equation 1711 complex eigenfrequencies, and therefore instability, are obtained, if V is finite and n is sufficiently 
large. For fixed n we obtain in the limit of large V (expansion of the root): 



2 
(7 2 



1 7r 2 n 2 
4 V 2 (a- b) 2 



2n 2 n 2 



V> 2 (a- bf 



(4-kt) 



(72) 



(73) 



Equation 1721 describes the eigenfrequencies of the decoupled shock modes discussed in section f4.1l i.e., the second order 
analysis of the previous section is contained in the limit V — > oo of the present approach. Instabilities described by equation 
I71l rcsemble those of the averaged model, rather than those given by equation [75] for positive values of k v . We conclude that a 
finite but large value of the stratification parameter V = d'ogp ^ s esseri tial for instability. However, assuming V — > oo, which 
was done in the investigation of the three-zone-model f section J4.1fl . is an oversimplification. 



5 CONCLUSIONS 

When following the non linear evolution of strange mode instabilities in the envelopes of massive stars, shock fronts were ob- 
served to be captured in the H-ionisation zone some pulsation periods after reaching t he non linear regime. This effect is not ob- 
served in models of very hot envelopes (such as the massive star model investigated bv lGlatzel. Kiriakidis. Chernigovskii fc FricI 
l)l99flFl 'l. due to hydrogen being ionised completely. The shocks trapped in the H-ionisation zone perform high frequency os- 
cillations (associated with the sound travel times across the shock zone) confined to its very vicinity, whereas the remaining 
parts of the envelope vary on the dynamical timescale of the primary, strange mode instability. By performing an appropriate 
linear stability analysis the high frequency oscillations were shown to be due to a physical instability, rather than being a 
numerical artifact. 

An analytical model for the secondary, shock zone instabilities has been constructed. As a result, high values of V 
were found to be responsible for instability. Contrary to the common stratification (convective, Rayleigh- Taylor) instabilities 
driven by buoyancy forces and thus associated with (non radial) gravity modes, however, the instabilities found here are 
associated with spherically symmetric acoustic waves. An extension of the stability analysis to non radial perturbations 
would be instructive, sin ce we expect the acoustic instabilities identified here - similar to strange mode instabilities (see 
Idatzel fc Mehreiil il996i) 'l - not to be restricted to spherical geometry. Such an investigation would also reveal buoyancy 
driven instabilities, which we believe not to be relevant for the following reasons: Their typical timescale is much longer 
than that of the acoustic instabilities, which will therefore dominate the dynamics. Moreover, in addition to gravity, the 
acceleration due to the shock's velocity field has to be taken into account and is likely to stabilize the stratification with 
respect to convective instabilities. With respect to the aim of this paper to identify the secondary, shock oscillations and 
their origin, a non radial analysis is beyond the scope of the present investigation and will be the subject of a forthcoming 
publication. 

Since the oscillations are a physical phenomenon - rather than a numerical artifact - they should not be damped by 
increasing the artificial viscosity as one would neglect a physical process whose influence on the long term behaviour of the 
system cannot yet be predicted. On the other hand, following the shock oscillations by numerical simulation for more than 
a few dynamical timescales is not feasible due to the small timesteps necessary to resolve them. The confinement of the 
oscillations to the very vicinity of hydrogen ionisation, however, indicates a solution of the problem by means of domain 
decomposition: The stellar envelope is decomposed into three domains: below, around and above the shock. Only the narrow 
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shock region needs high time resolution, the inner and outer zones merely require the dynamical timescale to be resolved. The 
development of a code following this strategy is in progress. 

Even if the appearance of the shock oscillations has so far prevented us from performing simulations in excess of several 
dynamical timescales, the velocity amplitudes reach a significant fraction of the escape velocity. This indicates that pulsation- 
ally driven mass loss may be found in appropriate simulations. Whether the new code will allow for the corresponding long 
term simulations and thus possibly for the determination of mass loss rates, remains to be seen. Preliminary results will be 
published in a forthcoming paper. 
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